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Abstract 

In  this  paper  we  discuss  the  sound  absorption  property  of  arrays  of  micro-acoustic  actuators  at  a  control 
surface.  We  use  the  wave  equation  over  the  half  plane  for  the  velocity  potential  with  a  boundary  dissipation 
by  a  proportional  pressure  feedback  law  along  the  half  plane  boundary.  The  feedback  gain  over  the  array 
is  described  by  a  distributed  shape  function.  We  develop  a  computational  method  based  on  the  Fourier 
transform  and  employ  it  for  analyzing  and  evaluating  the  decay  rate  of  acoustic  energy.  Specifically,  we 
carry  out  computations  for  a  diffusive  random  initial  field  and  report  on  our  resulting  numerical  findings. 


1  Introduction 


The  rise  and  decay  of  acoustic  energy  in  a  room  of  a  given  size  is  primarily  governed  by  the  area  and  absorption 
coefficient  of  the  surfaces  of  the  room  and  objects  within  [7,  10,  14],  The  absorption  coefficient  of  a  material 
with  normalized  input  impedance  (  =  6  +  i\  is 
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The  sound  absorbing  characteristics  of  passive  materials  are  excellent  when  the  wavelength  of  sound  is  less  than 
the  material  thickness,  because  the  input  impedance  is  largely  resistive.  However,  as  the  wavelength  of  sound 
increases,  the  reactive  part,  y,  of  the  impedance  increases,  and  the  performance  of  sound  absorbing  materials 
decreases  drastically  [6] .  The  apparent  solution  of  increasing  the  area  and  thickness  of  absorption  treatments  is 
often  limited  by  practical  considerations  including  cost. 

Active  sound  absorbers  provide  an  alternative  technique  for  achieving  sound  absorption  at  long  wavelengths 
where  the  size  of  the  absorber  is  on  the  order  of  the  wavelength  of  sound  [5,  8,  11].  Typical  active  sound 
absorbers  are  coupled  discrete  devices  or  systems,  for  example,  a  system  composed  of  a  microphone,  speaker, 
and  appropriate  control  law.  The  most  significant  limitation  of  active  sound  absorbers  is  their  limited  surface 
area,  and  again  practical  consideration  limit  their  use. 

Recently  the  concept  of  using  arrays  of  small  acoustic  sources,  which  could  be  constructed  of  fluidic  devices  [1] 
or  piezoelectric  materials,  for  sound  absorption  has  been  of  interest.  In  such  an  array,  each  element  would  be 
locally  reacting  by  feeding  back  a  pressure  measurement  via  an  appropriate  compensator  to  the  elemental 
acoustic  source.  The  dimensionality  of  an  array  of  practical  size,  however,  limits  the  analysis  of  the  array’s 
characteristics,  including  absorption,  using  techniques  used  in  previous  active  sound  absorber  analyses. 

In  this  investigation  we  develop  a  computational  technique  for  analyzing  the  absorption  rate  of  acoustic 
arrays.  Specifically,  we  consider  the  absorption  rate  of  an  ideal  reverberant  random  sound  field  as  defined  in 
Section  3. 

We  begin  with  a  brief  description  of  the  physical  system  under  consideration,  and  consider  a  typical  rate 
feedback  law  along  the  boundary  where  the  elemental  acoustic  source’s  volume  velocity  is  proportional  to  the 
measured  pressure  at  that  element.  This  corresponds  to  a  purely  resistive  input  impedance.  We  then  formulate 
the  acoustic  wave  equation  with  the  absorption  boundary  coefficient  in  the  frequency  domain.  In  Section  3  we 
describe  the  initial  sound  field.  In  Sections  4  and  5  we  formulate  the  system  in  the  frequency  domain  and  obtain 
a  discrete  model  approximation.  Based  on  this  discrete  model,  in  Sections  6  and  7  we  show  how  to  calculate 
the  instantaneous  correlation  field  and  energy  decay.  Finally,  we  report  on  a  numerical  experiment  and  analyze 
the  results. 


2 


y 

t 

Acoustic  Array 


0  - ►  * 


Figure  1:  Acoustic  array. 

2  Physical  Problem  Equations 

In  this  section  we  present  the  physical  system  under  consideration.  The  sound  pressure  p  satisfies  the  acoustic 
wave  equation 

Prr  ~  C2V2p  =  0,  (2) 

and  Euler’s  equation 

pvr  =  -Vp,  (3) 

where  c  is  the  speed  of  sound  in  air  and  p  is  the  mass  density  of  air.  The  operator  V  is  the  gradient,  v  =  ( u,v,w ) 

is  the  particle  velocity  vector  of  air,  r  denotes  time,  and  vr  =  dv/dr  ;  see  [4]  and  [13]. 

Consider  an  infinite  air  space  bounded  by  an  infinite  rigid  wall  located  at  x  =  0  covered  by  an  acoustic 
array  symmetric  about  the  origin.  We  assume  that  the  acoustic  array  is  long  enough  in  the  ^-direction  that  the 
system  is  independent  of  2  and  becomes  two  dimensional,  as  depicted  in  Figure  1.  Accordingly,  our  working 
space  is  ST  =  {(x,y)  :  x  >0,y  £  5ft}. 

The  interaction  of  a  wave  with  the  boundary  is  described  by  the  boundary  condition 

v(0,2/,t)  -i  =  u(0,y,T)  =  -g(y)p(0,y,T),  -oo  <  y  <  00,T  >  0.  (4) 

This  boundary  condition  follows  from  the  definition  of  the  normal  input  admittance,  the  ratio  of  the  particle 


3 


velocity  normal  to  the  wall  to  the  pressure  [6],  where  g  is  the  interaction  admittance.  Note  for  the  case  of  a 
rigid  boundary,  g  =  0,  while  for  a  pressure  release  boundary  we  have  g  — >  oo. 

Let  (f>  be  the  velocity  potential,  v  =  Vd>.  Then  from  (3)  we  have  p  =  —p<pT  and  it  follows  from  (2)  and  (4) 
that  0  satisfies 

<j>TT  -  C2V20  =  0  in  n,  (5) 

0x{O,  y,T)  =  p  g{y)  (j>r(0,y,T),  -oo  <  y  <  oo,  r  >  0.  (6) 

We  normalize  time  by  the  sound  speed,  letting  t  =  ct ,  and  we  define  7 (y)  =  pcg(y),  then  (5)  and  (6)  take  the 
form 


G 

0 

II 

> 

1 

(7) 

4>x(0,y,t)  =  7 (y)  4>t(Q,y,t), 

—  00  <  y  <  00,  t  >  0. 

(8) 

The  admittance  magnitude  in  a  given  direction  is  a  function  of  the  wave  propagation  angle.  Let  x  and 
k  denote  the  position  vector  (x,y)  and  the  wavenumber  vector  ( k,l ),  respectively.  Then  for  a  plane  wave, 
d>(x,  t)  =  <j)0  e'M-k-x)  5  the  normalized  wave  admittance  in  the  .T-direction  is 


ukk 
P  =  PC~  =  ~  =  Tj-r, 
p  u  k 


(9) 


since  cj  =  |k|.  We  note  that  0  <  |/?|  <  1  with  |/3|  =  0  corresponding  to  a  plane  wave  traveling  only  in 
the  (/-direction  and  |/i|  =  1  corresponding  to  a  plane  wave  traveling  only  in  the  .x-direction.  In  order  to  absorb 
sound  at  the  boundary  we  might  choose  the  normal  input  admittance  7  to  match  the  normal  wave  admittance  ft. 
However,  in  general,  the  knowledge  of  the  angle  of  incidence  of  the  incoming  wave  is  not  available.  Moreover,  for 
some  fields  the  angle  of  incidence  is  arbitrary.  For  this  reason,  in  our  model  we  test  the  absorption  characteristics 
for  a  given  7  by  the  absorption  rate  of  an  ideal  reverberant  sound  field  where  the  angle  of  incidence  is  uniformly 
distributed  among  all  directions. 


4 


3  Initial  Field 


The  sound  field  in  a  reverberation  room  is  assumed  to  have  the  property  that  it  is  random  and  at  every  point 
within  the  field  all  directions  of  the  incoming  sound  wave  are  uniformly  distributed  [3,  9,  12].  The  reverberation 
room  has  become  an  important  tool  in  applied  acoustics.  For  example,  it  is  useful  in  measuring  the  sound 
absorption  coefficients  of  surface  materials  and  the  sound  transmission  loss  of  building  structures.  Thus  in  our 
investigation  we  assume  the  initial  sound  field  to  be  random  and  the  spatial  autocorrelation  to  be  independent 
of  the  position  of  the  phase  and  a  function  of  the  distance  only.  That  is,  the  initial  random  field  0(x,  0)  satisfies 

£[<Kx,0)  0(x  +  r,O)]  =  i?(|r|),; 

for  all  x,  where  £  [•]  denotes  the  expectation. 

Let  <f>{ k,  0)  be  the  Fourier  transform  of  </>(x,  0)  , 

d>(k,0)  =  j  e_,k  x  d>(x,  0 )dx. 

Then,  as  in  [4],  the  frequency  autocorrelation  for  <^>(k,0)  is  given  by 

£  [<£(k,O)0(k',O)]  =  (2tt)2  6(k  -  k')  J  R(\r\)  e-*-r  dr,  (10) 

where  S  is  the  Dirac  delta  function. 

In  this  paper  we  develop  approximate  methods  to  evaluate  the  frequency  autocorrelation  for  the  sound  field, 

£  p(M)^(k ',*)]  , 

for  the  pressure  field, 

£  [<£t(k ,t)  &(k ',*)]  , 

and  for  the  velocity  field, 

£[v$(k ,t)  V^(k ',*)]  , 

as  well  as  the  total  energy  decay  of  the  sound  field  for  t  >  0. 
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4  Frequency  Domain  Equations 


In  this  section  we  reformulate  Equation  (7)  and  (8)  in  the  frequency  domain.  Since  the  spatial  domain  0  is  the 
half  plane  x  >  0,  then  without  loss  of  generality  we  use  the  Fourier  cosine  transform  with  respect  to  x  and  full 
Fourier  transform  with  respect  to  y.  Let  4>  be  the  Fourier  cosine-full  transform  of  <f>, 

/OO  cOO 

/  </>(x,  y,t)  cos(kx)  e~tly  dx  dy.  (11) 

-oo  J 0 

Then  the  inverse  transform  is  given  by 


1  poo  poo 

<f>(x.,  t)  =  <f>(x,  y,  t)  =  —  /  /  (j>(k,  l,  t)  cos(kx)  etly  dk  d.l.  (12) 

7T-  J_OQ  J o 

We  can  derive  an  equation  for  oik.  /,  /)  as  follows.  We  multiply  (7)  by  cos (kx)e~tly ,  and  apply  Green’s 
Formula  and  (8)  to  obtain 

/oo  coo 

/  <f>(x,  y,  t )  cos(kx)  e~tly  dx  dy  = 

-oo  J  0 

/OO  cOO  c  OO 

/  o(x.  y,  t)  cos (kx)  e~dy  dx  dy  -  /  'y(y)<t>t(0,  y,  t )  e~lly  dy.  (13) 

-oo«/0  J —oo 

Then,  by  applying  the  formulas  (11)  and  (12)  we  can  write  (13)  as 

/oo  poo 

/  < l>t(a,P,t)Af(l,P)dad/J ,  (14) 

-oo  J 0 


where 


1  P°° 

7(U)  =  ~  7  (y)e~,il 

^  J- oo 


-i(l-9)y 


Since  the  acoustic  array  is  assumed  to  be  symmetric,  7 (y)  is  an  even  function  and  thus 


7 (l,P)  =  4-  [  liy)  cos[(/  -  p)y\ 
7r“  J  0 


2  f 00 

dy  =  —  /  7(y)[cos(ly)  cos (0y)  +  sin (ly)  sin (py)]  dy.  (15) 


If  we  further  assume  that  the  initial  field  is  an  even  function  of  y,  then  the  map  y  (/>(-,  y,  •)  is  even,  and  we 
only  need  to  work  with  the  nonnegative  values  of  l.  Moreover,  the  integrand  7(2/)  sm(ly)  sm(/3y)  in  (15)  is  an 
odd  function  of  (3  and  thus  does  not  contribute  to  the  integral  term  of  (14).  As  a  result,  Equation  (14)  reduces 


4>tt(k,  l,  t)  =  ~(k’2  +  t2)(f>(k,  l,  t )  - 


4>t  ( a ,  Id ,  t )  T  ( l ,  /3)  da  did , 
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where 


4  r°° 

r  (l,P)  =  —  ~/(y)cos(ly)cos(py)dy.  (17) 

JO 

We  note  that  this  yields  an  integro  partial  differential  equation  for  (f>  for  which  approximation  techniques 
must  be  developed.  In  the  next  section  we  present  a  semi-discrete  finite  element  approximation  based  on 
piecewise  constant  elements  (zero-order  splines). 


5  Approximate  Model 


In  this  section  we  introduce  a  finite  dimensional  approximation  of  Equation  (16).  We  consider  the  truncated 


finite  domain 


{( k,l )  :  0  <  k  <  K,0  <  l  <  L}. 


(18) 


Let  0  =  ko  <k\  <  ko  <  ...  <  A:  v,  =  K  and  0  =  lo  <  li  <  h  <  ■  ■  ■  <  1&'—  L  be  partitions  along  the  k-  and  Z-axis 
in  the  ZcZ-plane.  Let  A  kj  =  kj  —  /q_i,  A  lj  =  lj  —  lj-i,  and  \.k;  x  (/,  \:lj  ■  Let  ( k-, .  I ;  i  denote  the 

midpoint  of  each  cell  Rij. 

We  apply  the  (piecewise  constant  in  frequency)  approximation 

M  N 

kk,i,t )  *  <i>”'  ik.i.i)  =  d9) 

i= 1  j=  1 


where 


\;r(k.n  = 


( k,l )  G  Rij, 


^  0  otherwise, 

and  $%jN(t)  is  an  approximation  of  /  ,.  /  ).  Then,  we  evaluate  (19)  at  the  midpoint  of  each  cell  Rjj.  This 
results  in  the  system  of  equations  of  the  form 


ij 


(t)  +  A 


<1.  : ' 


(t)  + 


(t)  =  0, 


(20) 


where 


Xjj  — 


kf  +  lj , 


9jn 


r  (ij,p)dp. 


In  —  1 


(21) 
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Note  that  we  should  use  the  superscript  MN  with  the  variables  in  (20)  since  they  depend  on  the  discretization 
in  (k,l).  However,  for  notational  convenience,  we  will  hereafter  suppress  this  superscript  whenever  no  confusion 
will  result. 

To  write  this  system  is  the  matrix  form  we  introduce  the  3iMA  vector 

$  =  [$11>  $21)  •  •  •  )  $Ml)  $12,  $22,  •  ■  ■  -  $a/2)  ■  ■  •  ,  $1jv,  $2JV;  ■  ■  •  ,  $mn]T, 


the  MN  x  MN  diagonal  matrix 

A  =  diag(AH,A^,...  '2  '2  '2  '2  '2  '2 

and  the  damping  matrix 


G  = 


We  define  the  M  x  M  matrix 


A*.  — 


\  2 

rl>  a12, 

X2 

a22 , ■ ' 

X2 

■  ■  ,  am2  i  ■ 

511 

512 

5ijv 

521 

522 

•  •  •  52.x 

5.x  i 

5jv2 

'  '  '  5  V  V 

Afci 

CN 

<1 

•  •  •  A  k 

Aki 

<1 

■■■  A  k 

Aki 

<1 

■■■  A  k 

(22) 


(23) 


and  let  D  be  the  Kronecker  tensor  product 


D  =  kron(G,Afc). 

Then,  the  second  order  system  (20)  can  thus  be  written  in  the  matrix  form 

$(i)  +D${t)  +  A$(t)  =  0.  (24) 


As  usual,  we  can  write  (24)  as  a  first  order  system.  To  that  end,  define  the  state  vector 


Then  (24)  takes  the  form 


Z  =  AZ,  (25) 

where  A  is  the  (2 MN)  x  (2MAT)  matrix  defined  by 

0  I 
-A  -D 

and  I  is  the  {MN)  x  {MN)  identity  matrix.  It  is  this  approximate  system  that  we  use  in  our  calculations 
reported  on  below. 

6  Instantaneous  Correlation  Field 

According  to  the  discretization  described  in  Section  5,  it  follows  from  (10)  that  a  natural  condition  to  require 
on  the  approximating  random  Fourier  components  {Z;(  0)}  is 

£[Zi(0)Zj(0)]=0,  if  *  ^  j.  (26) 

Define  the  corresponding  approximate  correlation  matrix 

C{t)=£[Z{t)ZT{t)].  (27) 

Then  (26)  implies  that  (7(0)  is  diagonal.  Thus,  we  can  write 

2  MN 

C(0)  =  w-elv)  (eW)  ,  (28) 

V=l 

where  e ^  is  the  unit  vector  in  3i2MJV  and  W  =  (w„)  is  the  vector  of  the  initial  field  discrete  power  spectral 
densities. 

From  (25)  and  (27)  we  find  that  C(t)  satisfies  the  matrix  Lyapunov  differential  equation  system 

C{t)  =  AC{t)  +C{t)  AT.  (29) 

Note  that  C  is  a  2MN  x  2MN  matrix  and  thus  we  have  a  system  of  size  (2 MN)2. 

The  system  (29)  is  a  large  coupled  system  to  be  solved  in  its  entirety.  However  the  solution  of  (29)  is  given 
by 

2  MN 

C{t)=  Y,  wvF^{t)(F^{t))  ,  (30) 

V=l 
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where  the  vector  F("'( t)  for  each  v  satisfies 


F{v){t)=  AF^{t),  F{V\Q)  =  .  (31) 

This  can  be  shown  as  follows: 

2  MAT 

c  =  J2  Wv  {fv  (i?M)T  +  f{v)  (^)t} 

v=l 
2  MN 

=  £  wv  {. AF <">  (f<">)  ' r  +  F<*>  (^F^)  "} 

V=l 

2MN 

=  Y  w"  {AF(V]  +  f^  (fM)t4t) 

V=1 

=  AC  +  C  A*, 

The  superposition  formula  (30)  is  especially  efficient  when  we  consider  fields  with  compact  support  Fourier 
components. 


7  Instantaneous  Total  Energy 

We  define  the  total  energy  of  system  (24)  by 

EMN  ( t )  =  +  E*IN  (t) ,  (32) 

where  the  “kinetic”  energy  and  “potential”  energy  are  given,  respectively,  by 

=  \j  |vr”(x,i)|2  dx,  Era)  =  \  j  [o;^(x./|]2  f/x, 

and  <j>MN{x.,t )  is  the  inverse  Fourier  transform  of  4>MJV(k,i)  defined  in  (19).  By  Parseval’s  theorem,  we  have 

Ev(t)  =  l*T(t)  \m,  EP(t )  =  ^T(t)m. 

We  remark  that  EyIN  is  an  approximate  measure  of  the  kinetic  energy  density  associated  with  the  velocity  field 
while  the  term  F“f  is  an  approximate  measure  of  the  potential  energy  density  due  to  the  pressure. 

Now,  by  multiplying  (24)  by  $  we  can  show  that 

E(t)  = -$T(t)D$(t)  <  0.  (33) 
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We  thus  see  from  (33)  that  the  system  is  dissipative  for  positive  7. 


To  describe  the  relation  between  the  total  energy  E  and  the  correlation  matrix  C,  we  define  the  matrix 


t/A 

0 

0 

I 

Then  from  (27)  we  have 


and  thus 


(\/A$)(\/A$)T 

\/A$$T 

<I(\/a$)t 

4 

1 _ 

£  [ E(t )]  =  jj£  [$T(i)  £  $(i)  +  $T(f)  $(*)]  =  ^  trace [£  C(t)  £]. 


(34) 


Our  goal  is  to  quantify  the  decay  mechanism  of  the  total  energy  as  t  — »  00.  By  integrating  (31)  numerically 
we  can  evaluate  F(v\t)  and  determine  the  decay  rate  of  each  mode.  Moreover,  we  can  form  the  correlation 
matrix  C  and  determine  the  decay  rate  of  a  field  composed  of  a  group  of  modes. 


8  Numerical  Implementation  and  Results 

In  this  section  we  demonstrate  how  the  above  approximations  and  resulting  computational  model  can  be  used 
to  investigate  the  energy  decay  rate  of  a  given  bandwidth  of  velocity  pseudo-modes  defined  by  (31). 

Consider  the  periodic  acoustic  array  of  length  2 La  in  the  (/-direction  and  symmetric  about  the  origin  as 
depicted  in  Figure  2.  The  array  consists  of  elements  each  of  width  w.  Because  of  symmetry,  we  only  need  to 
consider  the  upper  half  of  the  array.  The  number  of  elements  in  the  half  array  is  given  by 

Ne  =  La/  w. 

For  n  =  1, . . .  Ne,  each  nth  element  occupies  the  interval  Xn  =  [(n  —  l)w,  nwj. 

In  each  element  we  take  the  function  7 (y)  to  be  a  symmetric  triangle  curve  of  height  h  as  shown  in  Figure  2. 
Thus  we  have 

I  7 n(y),  </€!„, 

7  (y)  =  < 

0  y  >  La, 


11 


Figure  2:  Acoustic  array  and  the  admittance  function  7 (y). 


where 


2h 


7 n(y)  = 


y  —  ‘2(n  —  l)h,  y  £  [( n  —  l)w,  (n  —  .5)w], 


—  =^-y  +  2nh,  y  £  [( n  —  .5)w,nw]. 

For  our  example  calculations  we  use  La  =  4  meters  and  element  width  w  =  0.01  meters.  For  the  height,  h,  we 

choose  the  value  2.  For  this  value,  the  average  value  of  7  in  each  element  is  1. 

Take  I\  =  L  =  10  and  uniform  partition  A ki  =  A lj  =  1.  This  implies  that  M  =  N  =  10  and  that  A  of  (25) 
is  of  size  200  x  200.  We  calculate  the  entries  of  the  matrix  G  by  using  the  midpoint  approximation 


9in  «  A l  F (ij,ln)  =  —  A l  7 (y)  cos (Ijy)  cos (jny)  dy, 


(35) 


and  then  a  quadrature  rule,  for  example,  the  Trapezoidal  rule.  Figure  3  depicts  the  eigenvalues  of  the  resulting 
matrix  A. 

As  an  example,  we  consider  the  bandwidth  3.5  <  |k|  <  4.5.  The  discrete  modes  that  belong  to  this  bandwidth 
are  shown  in  Figure  4.  The  set  of  these  modes  is 

T  =  {(kijj)  =  ( *  —  -5,  j  -  .5 )  :  (t,j)  =  (4,1),  (4, 2),  (3, 3),  (4, 3),  (1,4),  (2, 4),  (3, 4)}. 


We  assume  that  the  initial  field  consists  only  of  the  velocity  field.  Thus  for  the  modes  in  T  the  corresponding 
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Figure  3:  Eigenvalues  of  the  matrix  A. 


Figure  4:  The  set  T  of  the  initial  discrete  modes. 
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set  of  v  =  i  +  10(j  —  1)  values  is 


Af  =  {4, 14,  23, 24,  31, 32, 33}, 


and  for  the  coefficients  wu  we  have 


wv  =  0,  is  ^  Af, 


while  wv,  v  G  Af,  are  given  magnitudes  that  describe  the  initial  field. 

As  a  result,  the  initial  correlation  matrix  (28)  reduces  to 

CA'(0)  =  Y,w,e(v)  (e(^)T  , 

veM 

and  the  corresponding  instantaneous  correlation  matrix  is  given  by 

C*(t)  =  (FW(i))T> 

v€Af 

where  is  the  pseudo-mode  given  by  (31). 

To  find  C^(t)  we  integrate  numerically  equation  (31)  for  each  F^\  is  £  AT.  We  use  the  MATLAB  function 
ode23s  which  is  suitable  for  the  stiff  system  (31)  Once  CA  is  calculated,  then  from  (34),  the  corresponding 
expected  instantaneous  total  energy  is  given  by 

£  [E-^(t)]  =  ^  traee[£ C^(t)C\.  (36) 


For  simplicity  we  consider  a  uniform  initial  energy  distribution  of  magnitude  one  for  each  mode  is  G  Af.  This 
implies  that 

wv  =  "  =  A jj~ ,  is  G  A  , 

where 


i  =  ind(i/,  10),  j  =  (is  —  i)/10  +  1. 


(37) 


Here  ind(a,  b )  is  defined  as 


ind(a,  b) 


b 

a  —  [[a./fo]]  x  b 


where  [[•]]  is  the  greatest  integer  function. 


if  b  is  a  divisor  of  a 
otherwise 
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Figure  5:  Evolution  of  the  energy  distribution  matrix  E(t)  at  t  =  0,2,4,  6,  depicted  in  (a),  (b),  (c),  and  (d), 
respectively. 


Since  CA  (t)  has  200  x  200  entries,  it  is  not  convenient  to  consider  a  graph  that  shows  the  simultaneous 
evolution  of  all  the  entries.  However,  since  the  significant  interactions  take  place  along  the  diagonal  entries 
of  ,  it  will  be  sufficient  to  consider  these  entries  only.  To  visualize  the  evolution  of  the  correlation  matrix 
associated  with  velocity  field  we  consider  the  entries  {  (CC^(t)  C)  ..j^.  We  form  the  10  x  10  energy  distribution 
matrix  SA  ( t )  by  mapping  these  entries  into  the  kl- space.  In  Figure  5(a-d)  we  present  the  bar  graph  of  the  matrix 
SA  (t)  at  t  =  0,  2,4,  6,  respectively.  The  solid  curve  in  Figure  6  represents  the  normalized  energy 


10  log10 


( £  [E^W]  \ 

\£[E*m  J  ■ 


By  examining  the  eigenvectors  of  A,  we  find  that  each  one  has  a  leading  entry  of  maximum  magnitude.  This 
enables  us  to  associate  each  eigenmode  with  a  pseudo-mode  F and  use  the  eigenmodes  to  obtain  an  estimate 
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Figure  6:  Energy  curves  for  K  =  L  =  10,  pseudo-modes  (solid),  eigenmodes  (dashed). 

to  the  correlation  matrix  and  total  energy  without  the  need  to  integrate  (31).  We  illustrate  this  as  follows. 

Let  ij}v  be  the  eigenvector  of  A  associated  with  the  eigenvalue  /A ,  that  is, 

=  HVV- 

Then  the  solution  of 

Tv{t)  =A  J-M(0 )  =  i/)v,  (38) 

is  =  exp (p,"f)  ipv.  We  call  an  eigenmode.  Now  we  can  use  the  set  of  eigenmodes  {T(v\t)} to 

calculate  an  estimate  to  the  correlation  matrix  and  the  expected  total  energy.  The  dashed  curve  in  Figure  6 
represents  the  normalized  energy  calculated  using  the  eigenmodes 

From  Figure  5  we  observe  that  the  energy  exchange  with  sufficiently  high  frequencies  is  negligible.  This 
indicates  that  it  is  not  necessary  to  choose  very  large  values  of  K  and  L.  In  Figure  7  we  compare  the  total 
energy  values  for  I\  =  L  =  10  with  M  =  N  =  10  (solid)  and  I\  =  L  =  20  with  M  =  N  =  20  (dotted).  As  can 
be  seen,  the  differences  are  negligible. 

Finally,  note  that  the  step  function  on  each  cell  [k,  k  +  Ak]  x  [l,  l  +  A l]  in  the  frequency  domain  corresponds 
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Figure  7:  Comparison  of  Energy  curve  for  K  =  L  =  10  (solid)  and  K  =  L  =  20  (dotted). 


to 

2  r  oo  /*oo 

<t>ki{x,y )  =  —  /  Xki(a,P)  cos(ax)  cos(fiy)  dadfi 

7T  ./o  ./o 

/s-|-A/s  /‘/-(-A/ 

/  cos(ax)  cos(fiy)  da  d/3 


in  the  physical  domain,  where  sine  .t  =  .  Thus  if  smaller  Afc  and  A l  are  chosen,  the  magnitude  of  (f>ki  is 

dispersed  over  a  large  area  and  leads  to  the  smaller  scale  interaction  of  the  sound  field.  Therefore,  we  must 
choose  an  appropriate  cell  size  A k  and  A l  in  order  to  admit  a  desired  scale  interaction.  For  example,  if  the 
spatial  domain  is  bounded  by  rigid  walls  of  size  Lx  and  Ly,  then  AA;  =  -f-  and  Al  =  -f-. 

We  remark  that  questions  of  convergence  to  a  stable  system  as  N,  M  — >  oo  in  the  problems  discussed  here  are 
not  particularly  relevant.  More  precisely,  the  margin  of  stability  of  the  original  system  (5)-(6)  is  zero.  In  fact, 
it  is  observed  in  our  calculations  that  the  distribution  of  eigenvalues  in  Figure  3  is  shifted  toward  the  imaginary 
axis  by  every  doubling  of  N ,  M .  One  can  employ  our  algorithm  in  a  more  practical  setting  by  restricting  the 


AkAl  .  A k 
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computational  domain  (i.e.,  the  model)  to  a  finite  cavity  with  weak  boundary  damping.  This  yields  an  acoustic 
problem  with  a  fixed  margin  of  stability  for  which  questions  of  limiting  behavior  as  N,  M  — y  oo  are  relevant 
(e.g.,  see  [2]). 

9  Concluding  Remarks 

The  above  discussions  constitute  our  initial  contributions  toward  the  effective  use  of  an  array  of  microacoustic 
actuators  as  an  absorbing  surface  for  enclosed  acoustic  cavities.  We  formulated  an  approximate  time  domain 
initial-boundary  value  problem  for  a  distributed  parameter  system  for  the  acoustic  energy  due  to  propagation 
of  an  initial  random  field.  Modeling  the  cavity  as  a  half-plane,  we  derived  the  corresponding  frequency  domain 
model  and  the  resulting  expressions  for  autocorrelation  for  the  sound  field,  including  the  pressure  field,  the 
velocity  field,  and  the  total  energy. 

We  then  developed  semi-discrete  finite  element  approximations  based  on  piecewise  constant  elements  (zero 
order  splines)  for  these  expressions  and  illustrated  their  use  with  sample  calculations. 

We  believe  the  computational  methods  discussed  in  this  paper  can  provide  useful  tools  for  predicting  the 
performance  of  a  given  actuator  array  with  respect  to  a  given  band  of  frequencies.  This  information  can  in  turn 
be  used  in  optimal  design  of  such  arrays. 
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